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Abstract. The time dependent surface flux (t-SURFF) method is extended to single 
and double ionization of two electron systems. Fully differential double emission 
spectra by strong pulses at extreme UV and infrared wave length are calculated 
using simulation volumes that only accommodate the effective range of the atomic 
binding potential and the quiver radius of free electrons in the external field. For 
, a model system we find pronounced dependence of shake-up and non-sequential 

■ double ionization on phase and duration of the laser pulse. Extension to fully three- 

dimensional calculations is discussed. 
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1. Introduction 

Differential double photo-electron spectra and the corresponding ionic recoil momentum 
spectra testify of dynamical correlation between the electrons. By sweeping extreme 
ultra violet (XUV) photon energies from below to above the threshold for single-photon 
double ionization of the He atom one probes correlation in initial and final states. In 
wave length in the infrared (IR) range, momentum distributions of recoil ions provide 
evidence for the importance of re-collision processes [1], where first one electron is 
ionized, which subsequently is re-directed by the oscillating laser field into a collision 
with its parent ion causing excitation and possibly detachment of the second electron. 
The early observation of unexpectedly enhanced double ionization of helium in IR fields 
[2] is now generally ascribed to this mechanism. Experimental data on strong field IR 
photo-ionization is also available for many other atomic and molecular systems and it 
was even proposed to use re-collision electron spectra for the analysis of structure and 
dynamics of molecules [3] . 

For the XUV wave-length, theoretical and experimental questions have matured, 
even if still under debate (see, e.g., [4] for a recent contribution to the debate with 
ample references to theory and experiment). At longer wave length, the large body of 
experimental data (sec, e.g. [5, 6, 7, 8, 9] and references therein) and the somewhat 
smaller range of theoretical models largely based on classical or semi-classical methods 
(see, e.g., [10, 11, 12, 13] ) are all plagued by the almost complete lack of reliable 
theoretical verification, with the notable exception of a few very large scale simulations 
of two-electron systems in strong fields [14, 15, 16], where, however, only in Ref. [15] the 
full two-electron dynamics is treated for laser wavelength of 800 nm. There are several 
reasons for this striking absence of complete ab initio simulations of ionization of two- 
electron systems. Firstly, even single ionization is non-trivial to compute, if the external 
fields are non-perturbative. Roughly speaking, the effort for computing single ionization 
grows with the 4th power of the wave length A"^ due to the growth of peak momenta oc A, 
quiver amplitude of the electron motion in the field oc A^ and the growth of minimum 
pulse duration ~ A (see also the discussion in Ref. [17]). When the effect of the field 
is perturbative, the situation for single ionization relaxes somewhat, as one basically 
only needs to know the initial neutral state and the single-electron stationary scattering 
solutions in the energy range of interest. Although obtaining scattering solutions may 
be difficult, there is a well defined procedure and the whole technology of electron- 
ion scattering theory available to approach the problem. For double ionization, also 
in the perturbative regime the situation is more complex. The convenient partition 
into bound and singly ionized spectral eigenfunctions cannot be continued to above the 
double-ionization threshold: eigenfunctions above the double-ionization threshold will in 
general have, both, single- and double-ionization asymptotics. For distinguishing single 
from double ionization we therefore invariably need the solutions at large distances. 
Even without the need for asymptotic analysis, scattering with open double ionization 
channels is a challenging task. 
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Numerical simulations, in principle, can provide the full answer. The asymptotic 
analysis is usually done by propagating the wave function until after the end of the 
pulse and then analyzing its remote parts either in terms of momentum eigenfunctions, 

1. e. plane waves, or, when the tails of the ionic Coulomb potentials are not considered 
negligible, in terms of two-body Coulomb scattering solutions. Three-body scattering 
solutions, which would obviate a purely asymptotic analysis, unfortunately, are not 
accessible. The by far largest part of the computational effort in these simulations goes 
into following the solution to large distances until the pulse is over and analysis can 
begin. It may be irritating to think that a large effort is made to simulate dynamics 
that is known exactly in the case of finite range potentials or approximately at sufficient 
distances from all Coulomb centers: the free motion of one or two electrons in a laser 
field. 

In a preceding publication [17] we have shown how, by absorbing the wavefunction 
and recording of flux before absorption, single particle photo-electron spectra can be 
computed using simulation volumes that only contain the relevant range of the atomic 
potential and the electronic quiver motion in the field. In this so-called time-dependent 
surface fiux (t-SURFF) method, asymptotic information is accumulated during time 
propagation rather then drawn from the full wave function after the end of the pulse. 
The equivalence of both approaches was proven in [17] mathematically and by comparing 
numerical results with literature. Depending on the system parameters and accuracy 
requirements, absorption radii for typical strong field setups can be kept as small as 20 
Bohr radii and as few as 90 linear discretization coefficients for the radial motion can give 
better that 1% accurate spectra over the complete spectral range. For Coulomb systems, 
due to the long interaction, larger radii of ~ 100 Bohr and 200 to 300 discretization 
points may be needed. 

In the present paper, we extend this approach to two-electron systems and 
multi-channel single ionization as well as double ionization spectra. A numerical 
demonstration is provided using a two times one-dimensional model system. As first 
physical results, we find phase and pulse-duration dependence of shake-up and double- 
ionization spectra, phenomena that are consistent with the spatial asymmetry of very 
short laser pulses and the re-collision mechanism for shake-up and non-sequential double 
ionization. 

2. The t-SURFF method for two-electron systems 

We first briefly summarize t-SURFF for the single particle case. The basic requirement 
is that there is some radius Rc beyond which the Hamiltonian reduces to an asymptotic 
one with known solutions. Let H(t) denote the time dependent Hamiltonian of our 
system and assume that there exists an exactly solvable Hy{t) that at large distances 
agrees with H{t): 

H^{t) = H{t) for \r\ > and Vi. (1) 
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For a single particle in a laser field described in velocity gauge and a short range potential 
V{r) = for \ f\ > Rc, Hy{t) is the Hamiltonian for the free motion in the laser field 

H^{t) = l[-iV-Ait)r, (2) 

where A{t) = — j^_^£{t')dt' for an electric dipole field £{t). Here and throughout we 
use atomic units h — rrie — — 1 unless indicated otherwise {rrie and e denote electron 
mass and unit charge, respectively, the Bohr radius results as the atomic unit of length). 
The TDSE with H^{t) has the Volkov solutions 

X,-(r-',i) = (27r)-^/2e-^*(^^'*)e*^ (3) 

i.e. plane waves times with time-dependence by the well-known Volkov phase 

^k,t) = ^ /* dt'[k-A{t')r. (4) 

Scattering describes the behavior of the time-dependent wave function ^(r, T) at long 
times T and large distances |f| > Rc- We choose T and Rc large enough such that at T 
the pulse is over 

H,{t)^-^A iort>T (5) 
and the wave function has split into its bound and scattering parts 

*(f,r) = *b(r,r) + *,(f,r) (6) 

with the properties 

*b(r,r)f«0 for|f1>i?c (7) 
^s(r,T)«0 for|f1<i?c. (8) 
The approximate sign in (7) refers to the tails of any bound state function, which decay 
exponentially with increasing Rc. The approximate sign in (8) refers to the fact that 
electrons with very low energies P/2 ~ may not have passed Rc at time T. We only 
need to analyze ^^5(7^, T) in terms of asymptotic functions Xg(r, T). As ^s(r, T) vanishes 
inside the radius R^ we can multiply the full ^(r, T) by the function 

— * 

write the emission amplitude for photo-electron momentum k as 

b{k,T) = {x,;{TMRc)\^{T)). (10) 
For obtaining a time integral over a surface flux, we write (10) as an integral of the time 
derivative and use the fact on the support of 9{Rc) both, Xjfc(^)0 ^(^)^) evolve by 
the same Hamiltonian Hy{t). We obtain 

{Xk{TMRc)\^iT)) = t r dt{xkmHv{t),9{Rcm{t)) (11) 

^0 

The commutator vanishes everywhere except on |r| = Rc'. the asymptotic information 
is obtained by integrating the time dependent flux through a surface at flnite distance 
Rc- For the further discussion of the single-particle case and numerical examples with 
short-range and Coulomb potentials we refer to Ref. [17]. 
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Figure 1. Partitioning of coordinate space into bound (B), singly ionized {S,S) and 
doubly ionized D regions. Single and double photo electron spectra are obtained by 
integrating the flux across the boundaries between the regions. 



2.1. Single-ionization into ionic ground and excited state channels 

The simplest extension of the single electron method is for computing single ionization 
into ground and excited state ionic channels. The complication is that in presence of a 
strong field the ionic state can differ from the field free ionic state due to polarization 
and one must make sure to count flux passing the surface into the correct ionic channel. 
In this section we consider only the lowest ionic states that remain bound and do not 
contribute to double ionization. 

For decomposing coordinate space into bound and asymptotic regions we define on 
both coordinates fi and r2 projector functions 9i{fi, Re) and ^2('^2, Rc), respectively, as 
in the single particle case (9). Again picking sufficiently large times T and a sufficiently 
large surface radius Rc we can partition the wave function ^'(T) into its bound, singly 
ionized, and doubly ionized parts (see Figure 1) 

ii-e,)ii-02)^{T) (12) 
+ ^i(i - e2)^{T) + (1 - 0i)e2^{T) + e^e^^iT) 

:= ^„(T) + ^,(T) + ^^(T) + ^rf(T) (13) 

Here and in the following we suppress the arguments ri,r2 and Rc of Oi and 62- Note 
that the assignment of singly and doubly ionized character to the different regions is 
asymptotically exact: the error can be made arbitrarily small for any specific solution 
^{T) by choosing sufficiently large T and Rc. 
The single ionization Hamiltonian 

Hs{t) = H,{t) ® Hionit) (14) 

agrees with the exact Hamiltonian on the support of ^1(1—^2) (region S in Figure 1). The 
corresponding Hamiltonian on S is obtained by particle exchange. Channel solutions 
Xc,ik(^i> ^2, t) for the TDSE on S 

^4xc(ri, r2, t) = Hs{t)xc{ri: r2, t) (15) 
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have the form 

Xcki^i, f2, t) = (27r)-^/^e-^*We^^^- ® Mr2, t). (16) 
Here 0c(^2,^) solves the ionic TDSE 

i^(f>c{r2, t) = Hio„{t)(l)c{r2, t). (17) 

Rather than an initial we use a final condition that 0c(^2,^) ior t > T should be 
the desired final ion state (we need to solve the TDSE backward in time). The 
channel solution on the support S of the particle exchanged projector (1 — 61)62 is 

With the spectral amplitude 

b{k,c,T) = (x,-,e(7^)l^i(l - ^2)|*s(T)). (18) 

— * 

the probability density for finding at time T an electron with momentum k in ionic 
channel c is 

a{k,c,T)=2\b{k,c,T)\\ (19) 

The factor 2 arises from adding the two identical exchange symmetric contributions. 
For converting this integral into a time integral over a surface we make the simplifying 
assumption that the ionic solution never leaves the bound area 

0c(r2,t)^O for \f\>Rc and Vt. (20) 

This is the precise version of the assumption that the ionic states considered does not 
get further ionized. The approximate sign refers to the fact that again there is always 
an exponential tail reaching to arbitrary distances and further that interaction with any 
pulse will lead to a small amount of ionization. Using (20) we neglect the flux between 
the singly ionized regions S, S and the doubly ionized region D at all times. Then, using 
the same procedure as for a single particle we obtain 

b{k,c,T)^t [ dt{k,t\[H,{t),6,]\iP,{t)) (21) 
J —00 

where ipdrijt) is the channel projected wave function 

A{n,t):^ J d^r2(p:{f2,t)^{n,f2,t). (22) 

For evaluating the integral (21) we only need to know values and radial derivatives of 
■0c on the surface |r| = Re- 

For the computation it means solving the full two-electron problem up to time T 
and up to radius Rc- Beyond Rc one can absorb all amplitudes. In addition, for each 
channel c, we need to solve one single electron problem up to the same time and radius 
(which is usually a much simpler task). 



t-SURFF: two-electron photo emission 



7 



2.2. Double ionization spectra 

When double ionization occurs, flux passes from the bound region B through the singly 
ionized regions S or S into the doubly ionized region D. Our naming of the areas is 
suggestive but does not imply any bias as to an actual state that the electrons occupy 
within any of these areas. It is unsubstantial for the present discussion whether some 
intermediate ionic bound state is occupied by one of the electrons in S or S (sequential 
ionization) or whether both electrons must be considered unbound. The sole purpose 
of the partitioning is to have well-defined surfaces outside the ranges of the respective 
potentials where we will integrate fluxes. 

For double ionization there is one obvious limitation of the discussion so far: on 
the line |ri — = a the electron-electron interaction is constant and not negligible 
for small a. This problem is not related to the long-range Coulomb potential, rather it 
is a general problem for any multi-particle breakup, which is why break-up processes 
are more complex than single particle scattering. Within the framework of the present 
approach the problem can be solved for short-range electron-electron repulsion without 
making approximations, which will be discussed below. A pragmatic solution has been 
effectively employed in many earlier publications, which is to neglect electron-electron 
repulsion at large distances from the nucleus. This is what using any projection onto 
products of single electron states implies, be it Coulomb scattering waves or plane waves 
(both approaches are discussed, for example, in [18]). Sensitivity to this approximation 
can be tested by varying the distance from the nucleus where the projection starts. 

As the pragmatic solution was found to work well in many cases, we make this 
approximation explicit in the present paper by smoothly turning off all potentials 
including the electron-electron interaction before the surface radius Re- In that case 
we can always use the free (Volkov) Hamiltonian H^(t) beyond R^. 

By our assumptions, in the region |fi| > Rc, the Hamiltonian is identical to Hs{t), 
Eq. (14), which motivates the ansatz 

^i*(ri,f2,t)- [ d'kJ2Xki^i,t)Ur2)b\k,n), (23) 

with the Volkov solutions X^(ri,t) on coordinate fi and an expansion into a time- 
independent, complete, but otherwise arbitrary set of functions Cni'^i) on rg. Using 

— * 

orthogonal projection onto the expansion functions ^-nd the coefficients b'{ki, n, t) 
are obtained as 

b'{ki,n,t)^ J d^k[qe{h,k[,t)b{k[,n,t) (24) 

with 

/oo /»oo 
d'nxl,{fi:t)e, / C(^"'2)*(ri,f2,i). (25) 
■oo ^ J —oo 

The integral (24) over qo accounts for the fact that the plane waves are not (5-orthonormal 
when the integration is restricted by 9i: the inverse overlap is defined by 

J d^k'q0{k,k',t){k',t\ei\k",t) = 5^^\k-k"). (26) 
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For notational brevity, we assume here that the basis functions ^„ are orthonormal. The 
time-derivative of the b{ki,n,t) is 

ij^b{h,n,t) = {x^,^{t)\e,{umt)\^m 

-{x-,^it)\mt),oi]{U'^m (27) 

By the double right bracket )) we emphasize that integration is over both coordinates fi 
and r2. Inserting the representation (23) into the first term we obtain an inhomogeneous 
equation for the b{ki,n,t): 

i^b{ki,n,t) = '^{^n\Hion{t)\m)b{ki,m,t) 

m 

-{h,t\[H,{t)M{in\'^m (28) 

where we have used (26). The inhomogeneity is the flux through the surface = R^.. 
Initial conditions are b{ki,n) = 0, i.e. no electrons outside Re- We write the double 
ionization amplitude 

bik„k„T) = {XR,{T)\{x^^{T)\e,e,mm (29) 
and use one more time the conversion to integrals over surface flux 

b{h,h,T) = (x,^^(T)|(x,^^(T)|Mi|M/(T))) (30) 

= f dtj^{x^,^mxuM^2eA^i^))) (31) 

=: i J dt[B{ki, k2, t) + B{ki, h, t)]. (32) 

The two terms B and B are related by exchange symmetry 

B{h,k2,t) = B{k,,h,t) = {x^,^mXkMHvit),e2]e,\<iim (33) 

For computing B, we only need to know ^i^(t), for which we insert the representation 
(23) to obtain 

B{h,k2,t) = Y.(Xk,itmv{t),92]\UKk2,m,t) (34) 

m 

The inverse overlap qe{ki,k[), (26), cancels with the overlap integrals and never needs 
to be evaluated explicitly. 

B{ki,k2,t) is the contribution to the double ionization spectrum passing at time 
t through surface |r2| = Rc from region S into D. In S, the first electron is already 

— * 

detached and has a fixed canonical momentum ki that is carried into the double-ionized 

— * — * 

region D. B{ki, k2, t) is the alternate contribution going through region S. 
2.3. Computational remarks 

The substantial gain of the method is that, rather than computing the full solution in 
region D and then analyzing it, we only need to integrate the flux through the surface 
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separating S from D. In the direction parallel to that surface the wave function is 
represented in terms of the free solutions Xfe^l^i) where we do not need to expand the 
wave function completely, but can restrict propagation to the momenta ki that we are 
interested in. For each ki we need to solve equation (28), which is an ionic TDSE with 
an additional source term (the flux entering S from the bound region B). Although the 
derivation may appear complex, the implementation of the procedure can be done by 
the following simple algorithm 

• Set up b{ki,n,t) on a ki grid of the desired density and initialize to 0. 

• Get time grid points with a time step that resolves the fastest amplitude 
oscillations for the desired energy range At < 271 / E max 

• In a loop through all times tj, get the surface terms from file or from a simultaneous 
calculation of |^(tj))) and use (28) to advance to b{ki,n,ti). 

• Add the contribution from region S to the spectral amplitude 

— * — * — # — # — * — * 

s{ki, k2, ti) = s(/ci, A;2, U-i) + AtB{ki, k2, t). (35) 

The contribution s from region S is obtained from s by particle exchange. The total 
spectral amplitude is the sum of both contributions 

6(^1, T) = s(^i, k2, T) + s(^i, 4, T) = s(^i, ^2, T) + s(4, k^, T) (36) 

The asymptotic value — the spectral amplitude — is attained at times T when the flux 
through the surfaces becomes neghgible. 

The approach can only be successful, when absorption does not significantly distort 

the solution at the integration surfaces. The "infinite range exterior complex scaling" 
(irECS) absorber was shown in Refs. [17, 19] to provide traceless absorption over a 
very wide energy range at low computational cost. It outperforms standard complex 
absorbing potentials by several orders of magnitude in accuracy. If needed, irECS can 
be pushed to full machine precision using not more than 20 discretization coefficients 
per coordinate in the absorbing region \ f\ > Aq. Absorption can begin at any Aq > Rc. 
More mathematical and numerical detail and ample numerical examples using irECS 
can be found in Refs. [17, 19]. 

In general discretization errors for two electrons are similar to those for a single- 
electron system with the same ionization potential. This has the practical advantage 
that a suitable discretization can be determined from the single-electron problem. Only 
a few consistency checks on specific two-electron observables need to be performed for 
the computationally heavier two-electron calculation. 

2.3.1. General single ionization spectra Once we have obtained the b{ki,n,T) we can 
reconstruct the wave function for \fi\ > Rc. In particular, the amplitude for single 
ionization spectra for any ionic state 0c is 

b{k,c) = J2KKn)al^\ (37) 

n 
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where On"* are the expansion coefficients of 0c with respect to the basis functions 

00(0 = 5^^(04^). (38) 

n 

This approach is not hmited to non-ionizing (pc, as it analyzes ionic population at time T 
after the end of the pulse. For non-ionizing 0^, it is an alternative to the single-ionization 
procedure above. Note that here we need to solve the inhomogeneous ionic problem (28) 
for each photo-electron momentum k. The total spectrum is a linear combination of the 
individual contributions from each n. Where it is applicable, the advantage of the single- 
ionization procedure of section 2.1 twofold: firstly, for each final (pc one needs to solve 
only one ionic TDSE and compute values and derivatives of the channel surface function 
ipc- The complete spectrum can be obtained by time- integration with little numerical 
effort. Secondly, as one directly obtains the channel spectrum without intermediated 
decomposition and final resummation of the wave function, results are in general more 
robust numerically. 



3. Numerical demonstration of the method 
3.1. The one- dimensional two- electron model system 

As even with the simplifications by t-SURFF the solution of the full three-dimensional 
two-electron problem remains a very large scale computational task, we use for 
demonstration purposes the standard one-dimensional two-electron model Hamiltonian 

For making all potentials strictly finite range, wc have chosen a "truncation radius" Cp 
and a truncation function M{x) that is = 1 up to \xa\ = Cp — 5 and goes differentiably 
smoothly to at Ixq,] = Cp. Where not indicated otherwise we use Cp — 20. The 
screening factor of 1/2 in the ionic potential was chosen for esthetic reasons, as then the 
exact ionic ground state (without potential truncation) is Eq — —2. The first excited 
state occurs at Ei — —0.93, substantially stronger bound than the first excited He'^ 
energy of —1/2. With the electron-electron screening of 0.3 one obtains an ionization 
potential of 0.88, which we consider a fair approximation of the actual He ionization 
potential of 0.90. 

For all computations we use a finite elements discretization where orders between 8 
and 20 where used for convergence studies. For the present purposes we found order 8 
or 10 sufficient. A deeper discussion of the finite element method for irECS calculations 
can be found in [19]. The irECS radius Aq which marks the beginning of absorption 
was varied between 20 and 100 atomic units. Results do not depend on Aq on the level 
of accuracy shown. 



2M{xo 



+ 



M{xi)M{x2] 



y/xl + 1/2 A/fxi - xo)^ + 0.3 



.(39) 



t-SURFF: two-electron photo emission 



11 



In all calculations below we use cos^ pulse shapes defined in terms of the vector 
potential 

S Tit 

A{t) = cos2(— ) siniut + (t>cEo), (40) 

U Zl FWHM 

where TpyyHM is the full- width-half- maximum of the envelope of the vector potential, u 
is the laser central photon energy, and (f)c'EO the carrier-envelope offset phase. The peak 
field amplitude Sq is related to the pulse peak intensity / by = \/27. By construction, 
the electric field S{t) = has no zero frequency component, which is important 

when studying effects of 0cbo with very short pulses. For a "cosine pulse" 0c£O = 
the peak of the electric field approximately coincides with the maximum of the pulse 
envelope with minor deviations due to the derivative of the envelope. 

3.2. Two photon double ionization in the extreme ultraviolet 

In Ref. [4] perturbative two-photon double ionization of He by short XUV pulses 
with photon energies between the two- and single-photon ionization thresholds was 
investigated and a remarkable universal description of the process was found. The 
findings were confirmed by numerical solutions of the fully three-dimensional two- 
electron TDSE. Pulse durations of T = 4.5 fs and photon energies u between 42 and 
80 eV were used. The perturbative regime at = 42 eV extends up to intensities 
near / = lO^^W/cm'^, where one first observes sizable deviations from the perturbative 
scaling oc / at the single-photon single-ionization peak and oc P for shake-up and double 
ionization. 

As an example. Figure 2 shows the two-electron spectrum for photon energy uj — 
70eV at peak intensity lO^^W/cni^. Three two-electron at energies |£'i| + |i?2| = uuj + Eq 
are clearly visible. With ground state energy Eq = — 2.88 au ~ —78eV the first peak 
appears for n=2-photon ionization. The ridges parallel to the energy axes are artefacts 
of the method: there is significant shake-up and the exponential tail of the excited ionic 
states reaches into region D. Note, however, that with the choice of Rc = 40 for the 
present calculation, the ridges are suppressed by at least a factor 10~^ relative to the 
surrounding signal. Even stronger suppression can be obtained by further increasing Rc 
or, alternatively and with slower convergence rate, propagating to longer times. If the 
ionic states are known accurately, their effect can be completely removed by projection 
(for a detailed discussion of the procedure see Ref. [17]). 

In Figure 3 we show line-outs of the two-photon two-electron energy ridges for 
photon energies between 42 and 80 eV. The curves are converged within the resolution 
of the plots with respect to T, Rc and spatial discretization. Although the general 
characteristics are determined by the single-excitation resonances as in [4], for the 
present system we do not reproduce the universal behavior found there for the angle- 
integrated spectra. The likely reason are shake-up interferences that also modify the 
the angular distributions in the 3D case [4]. 

Even at the present benign, nearly perturbative parameters, photo-electrons with 
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Figure 2. Double electron emission by a 4.5 is XUV pulse at photon energy u = 70eV. 
with peak intensity IQ^^W/cm?. Negative energies indicate emission to (— oo,— i?c]- 
Straight lines at 45 degrees with constant Ei ± E2 indicate correlated processes. The 
first 3 double-electron ridges are visible. Lines parallel to the energy axes are artefacts 
due to exponential tails of excited ionic bound states that reach beyond Rc = 40. Their 
peak values lie below < 10~^ of the peak signal. Black color (< 10"^^°) approximately 
indicates numerical noise. 




Figure 3. Two-photon double ionization for photon energies 42, 54, 58, 70, and 
80 eV. Lineouts are taken at the respective lowest energy ridges with total energies 
l^-il + |-E'2| = 6,18,30,38,62 and 82 eV^, respectively. Curves are normalized to 1 at 
El- E2 = 0. 
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energies of about ~ 3au ~ 80 eV travel to a distance r^a-^ ~ 2TFWHM\/'^Ei^/me ~ 
900 aw during the duration of the pulse. To correctly represent the corresponding 
momenta one needs a grid spacing of Ar < 2ti / \/2Emax ~ 2.5, which results in at 
least about 400 radial discretization points in a numerical calculation. In practice, finer 
grid spacings need to be used with a correspondingly larger number of discretization 
points. In addition, for asymptotic analysis of the wave function, one may need to 
propagate further for a certain period of time after the pulse, leading to further increase 
in the required box radius Vmax proportional to propagation time. 

For comparison, with t-SURFF we obtain converged results using only 49 
discretization coefficients on the positive half axis [0,oo). The total number of 
discretization points per coordinate is twice as large because of the negative half- 
axis, which would correspond to a different angular direction in a three-dimensional 
calculation. 32 out of the 49 points were used on the interval [0,20], i.e. an effective 
grid spacing of ~ 20/32 = 0.625, which sets a theoretical limit for the maximum photo- 
electron energy of < 50aii fa 1300 eF. The remaining 17 points were used for irECS 
absorption. The t-SURFF simulation volume radius Rc is independent of pulse duration. 

3.3. Single ionization and shake-up by an IE pulse 

Possibly the most elementary correlated process is shake-up: after forcefully removing a 
single electron, the remaining electrons rearrange not exclusively into the ionic ground 
state, but a fraction goes into excited states. The exact distribution of excited states 
strongly depends on system parameters. Shake-up by very strong, short infrared pulses 
has been observed recently using the "attosecond transient absorption spectroscopy" 
[20], which can monitor the time-evolution of an excited ionic state after IR ionization. 

Within the present simplified model wc can only demonstrate basic qualitative 
features of shake-up in IR ionization. Wc have studied dependence on the carrier- 
envelope offset phase (I)ceo and on pulse duration. We use an intensity of 2 x lO^^W/cm^ 
and uj — 0.057au (corresponding to wave length 800 nm). Converged results were 
obtained using the same discretization parameters as in the XUV case discussed above. 

In Figure 4 photo electron spectra for the ground and first excited states calculated 
by formula (21, smooth line) and, alternatively, by Eq. (37, coarser energy grid) for a 
single cycle cosine pulse: T — 2t:/ui and (pcEO = 0. The spectra show some generic 
short pulse features: pronounced asymmetry and absence of individual photo-electron 
peaks at energies, where emission occurs only during a single laser half-cycle. Emission 
to the left is lower, as the field amplitude to the left is smaller. 

Clearly, the spectral structure is highly sensitive to (pcEO- However, also the total 
shake-up yield is strongly (pcEO dependent. Table 1 list the yields in the ground and 
first excited ionic states for different phases and pulse durations. While shake-up is 
strongly suppressed for a single cycle cosine pulse, it increases to a sizeable ~ 10% of 
the ground state ionic channel for (pcEO = 37r/4. We tentatively ascribe this fact to a 
recoUision mechanism for shake-up. The fact that longer pulses with 2 and 3 optical 
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Figure 4. Single photo-electron emission spectra in the lowest 4 ionic channels. 
Negative energies indicate emission to (— oo,— i?c]- Spectra are calculated using 
Eq. (37). In addition, for ground and first excited ion channel, curves computed by 
Eq. (21) arc shown (black lines). The curves are scaled for better visibility. 



TpWHM 


4>CEO 






yi 








1 opt. eye. 





2.26 X 10~ 


-7 


6.27 X 10- 


10 


2.77 X 10" 


-3 




7r/4 


2.18 X 10- 


-7 


1.35 X 10- 


-8 


6.19 X 10- 


-2 




7r/2 


8.88 X 10- 


-8 


1.35 X 10- 


-8 


1.15 X 10- 


-1 




37r/4 


1.38 X 10- 


-7 


2.60 X 10" 


-9 


1.88 X 10- 


-2 


2 opt. eye. 





2.36 X 10- 


-7 


1.16 X 10- 


~8 


4.91 X 10- 


-2 


3 opt. eye. 





1.44 X 10- 


-7 


3.18 X 10- 


-9 


2.21 X 10- 


-2 



Table 1. Shake-up as a function of 4>ceo- Total yield in the ionic ground Yq and first 

excited Yi states and ratio Fq/ bi- 
cycles in general have shake-up fractions comparable to the (t)cEO = 37r/4 is consistent 
with this hypothesis. We abstain from a more detailed analysis of the phenomenon for 
the present model, as due to the absence of transverse wave-packet spreading, recollision 
effects are greatly exaggerated in one dimensional models. 

3.4- Double ionization spectra generated by an IR pulse 

Figure 5 gives an overview of doubly-differential photo emission spectra obtained for 
different (pcEO and pulse durations up to three optical cycles Tfwhm = Qn/ojau ~ 
7.5 fs. The calculations were performed using a slightly more accurate discretization 
with 60 points on the half-axis [0, oo). 

As to be expected, we see pronounced asymmetries and strong effects of (pcEO ior 
the shortest pulses. While the single-cycle cosine pulse essentially only shows weak and 
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Figure 5. Overview of double ionization spectra for laser wave-length 800 nm and 
peak intensity 2 x . Upper row pulse duration T = 1 optical cycle with 

4'CEO — (left) and 3tt/4 (right). Lower row (peso = with T = 2 (left) and 3 
(right). Negative energies indicate emission to (— oo, —Rc\- 

unidirectional emission of both electrons, backward emission and significantly higher 
yields arise at (pcEO = 37r/4. Most likely, the responsible mechanism is re-collision. 
However, due to the limitations of the one-dimensional model, we defer a more profound 
discussion of the effects to a future fully three-dimensional calculation. 

3.4-1 ■ Dependence on the potential cutoff By truncating all potentials at Cp = 20 au, 
the Volkov Hamiltonian H^{t) is exact outside Rc and all errors of t-SURFF are either 
due to discretization or absorption. As a pragmatic approach such a truncation is 
very appealing, but the truncation error must be controlled. For our simple model, we 
made a series of calculations for truncation radii up to Cp = 100. Figure 6 shows the 
single photo-emission spectrum for the first excited ionic channel, a diagonal lineout 
El = E2 of the double-emission spectrum, and two lineouts where one energy energy 
is fixed at Ei = 0.1 0^(2.761^) and Ei = 0.5 a'u(13.6 eV^). The single electron spectra 
are hardly affected by the truncation. The strongest effect, as to be expected, is along 
the line Ei = E2. With Cp = 20 qualitative differences in the two-electron plots can be 
observed, while good qualitative agreement can be found for Cp = 50. The values Cp 
where results are acceptable must be chosen depending on the system and on accuracy 
requirements. In our case, the on-diagonal spectra are still incorrect up to Cp = 50. 
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Figure 6. Dependence of spectra on the potential cutoff Cp = 20, 50, 80 and 100 for a 
800 nm single-cycle cosine pulse with peak intensity 2 x IQ^^W/cm? . Top left: single 
ionization in the first shake-up channel ai{E)^ top right: double ionization with equal 
energies cf{E^ E), bottom: lineouts of the double ionization signal (j{E, Ei ) for energies 
Ei=3eV (left) and 13 eV (right). 
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Figure 7. Integrating the flux for non- vanishing e-e interactions. 



In contrast, convergence for larger momentum differences appears quite acceptable at 
Cp ^ 50. 

Alternatively to increasing Cp, one may use a different arrangement of surfaces, 
where an additional surface is placed at \ fi + = Rc (see Figure 7). In scaled center- 
of-mass and interparticle coordinates r± := (n + r2)/V^ the asymptotic Hamiltonian 
H£,{t) in region D separates into the center-of-mass motion with a Volkov solution and 
field-free relative motion in the repulsive potential of the electrons 

Hoit) = \ [-zV+ - V2i(t)l' - + ree(V2|r-|), (41) 



2 



2 



where V+ and A_ denote the Nabla and Laplace operators in coordinates Tj^ and f^_, 
respectively. Thus the asymptotic problem reduces to field-free potential scattering and 
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free motion in the laser field. The contribution going from region B directly into D can 
be conveniently expanded in this form. We leave the technical discussion and numerical 
demonstration of this procedure for future work. 

4. Discussion and conclusions 

With the reduction of the simulation volume afforded by t-SURFF, fully dimensional 
calculations of double photo-ionization for a broad range of wave-length and intensities 
has come into reach. Scaling to three dimensions depends on the laser frequency: in the 
perturbative regime, the number of angular momenta required may remain as low as 3 
[18], which, in linear polarization, scales the problem size by ~ 3^ = 27. In that case 
the computational effort is less than one order of magnitude larger than the calculations 
presented in section 3.2, if one considers that the positive and negative half-axes of the 
one-dimensional model lead to a 4-fold scahng of the problem compared to a radial 
problem. At infrared wave length, the complete angular phase-space becomes activated. 
At 800 nm wave length and intensities of ~ 2 x 10^^ VF/cm^ 30 angular momenta are 
needed for full convergence of electron spectra. For two-electron spectra, this entails 
an increase in problem size by ~ 30^/4 7000 compared to the present calculation, 
if we keep the same standards of accuracy. Considering that a typical solution for our 
problems takes 2 hours on a single CPU, one sees that 3D calculation are quiet feasible 
on a moderate size parallel computer. The fact that in the 3D case the fraction of the 
phase space, where electron repulsion is important, is much smaller than in 1-d, may 
possibly be exploited for further reduction of the problem size. 

The development of t-SURFF is motivated by the photo-ionization problem. What 
is special about that problem is that in dipole approximation the time-dependent field 
affects the whole wave function, including the asymptotic region and therefore final 
momenta cannot be determined before the end of the pulse, unless one uses knowledge 
about the time-dependence of the asymptotic solution. However, also in situations 
without time-dependence of the asymptotic Hamiltonian, the method may turn out to 
be useful for obtaining fully differential momentum spectra. All it requires is a reliable 
absorption method and knowledge of the solution in the asymptotic region. Among 
the candidates for further application are reactive scattering and chemical break-up 
processes. 
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